#include<iostream>
#include<vector>
#include<queue>
#include<malloc.h>
#include <math.h>

void FMG_cycle_one(double *f,double RA,double *u,int ng,int n,int RO,double *au,int IO,int ncycle_M,double *y)

{
  void addint_one_l(double *uf, double *uc, double *res, int nf);
void addint_one_q(double *uf, double *uc, double *res, int nf);
void copy_one(double *aout, double *ain, int n);
void fill0_one(double *u, int n);
void interp_one_l(double *uf, double *uc, int nf);
void interp_one_q(double *uf, double *uc, int nf);
void relax_one(double *u, double *rhs, int n);
void resid_one(double *res, double *u, double *rhs, int n);
 void rstrct_one_fw(double *uc, double *uf, int nc);
void rstrct_one_i(double *uc, double *uf, int nc);
void slvsml_one(double *u, double *rhs);
 double maxnorm_one(double a,double *u,int n);

 int nn,jj,ngrid,j,nf,jpost,icycle,jpre,i;
 double  *ires[16],*irhs[16],*iu[16],*irho[16],*mu,*er;
 double Maxnorm,umaxnorm,a=-1,h,h2i,pumaxnorm,rateR;
 
 h=1.0/(n-1);
 h2i=1.0/(h*h);
  mu=(double *)malloc(sizeof(double)*2000);
  er=(double *)malloc(sizeof(double)*2000);
for (icycle=1;icycle<=ncycle_M;icycle++) {
  std::cout<<"The "<<icycle<<"-th FMG cycle:"<<std::endl;
  std::cout<<" "<<std::endl;
   
  
 //ires,irhs,iu,irho全部给空间
  for (jj=ng;jj>=1;jj--) {
    nn=pow(2.0,jj)+1;
    irhs[jj]=(double *)malloc(sizeof(double)*2000);
    ires[jj]=(double *)malloc(sizeof(double)*2000);//残差
    iu[jj]=(double *)malloc(sizeof(double)*2000);
    irho[jj]=(double *)malloc(sizeof(double)*2000);
    fill0_one(iu[jj],nn);fill0_one(ires[jj],nn);fill0_one(irhs[jj],nn);
  }

  nn=(n+1)/2;
  ngrid=ng-1;
 
  if (RO==0)
      rstrct_one_fw(irho[ngrid],u,nn);
    else
      rstrct_one_i(irho[ngrid],u,nn);
  while (nn>3){
    nn=(nn+1)/2;
    ngrid=ngrid-1;
   
     if (RO==0)
      rstrct_one_fw(irho[ngrid],irho[ngrid+1],nn);
    else
      rstrct_one_i(irho[ngrid],irho[ngrid+1],nn);
  }
  slvsml_one(iu[1],irho[1]);
 
  ngrid=ng;
 
  for (j=2;j<=ngrid;j++){
    nn=pow(2.0,j)+1;
    if (IO==0)
      interp_one_l(iu[j],iu[j-1],nn);
    else
      interp_one_q(iu[j],iu[j-1],nn);
    copy_one(irhs[j],(j !=ngrid ? irho[j] :u),nn);
	nf=nn;
	//V_cycle begin
	for (jj=j;jj>=2;jj--){
	  for (jpre=1;jpre<=3;jpre++)
	    relax_one(iu[jj],irhs[jj],nf);
	  resid_one(ires[jj],iu[jj],irhs[jj],nf);
	  nf=(nf+1)/2;
	  if (RO==0)
            rstrct_one_fw(irhs[jj-1],ires[jj],nf);
          else
            rstrct_one_i(irhs[jj-1],ires[jj],nf);
	  fill0_one(iu[jj-1],nf);
	}
	slvsml_one(iu[1],irhs[1]);
	nf=3;
	for(jj=2;jj<=j;jj++){
	  nf=2*nf-1;
	  if (IO==0)
            addint_one_l(iu[jj],iu[jj-1],ires[jj],nf);
          else
            addint_one_q(iu[jj],iu[jj-1],ires[jj],nf);
	  for(jpost=1;jpost<=3;jpost++)
	    relax_one(iu[jj],irhs[jj],nf);
	}
	//V_cycle end
  }
  if (icycle==1)
   pumaxnorm=maxnorm_one(a,u,n);
  else
    pumaxnorm=umaxnorm;
  //重新计算精确值和残差
 
  for (i=2;i<n;i++){
    au[i]=au[i]+iu[ng][i];
  }
  for (i=2;i<n;i++){
   u[i]=-h2i*(au[i-1]+au[i+1]-2*au[i])+f[i];
  }
   u[1]=-h2i*(au[2]-2*au[1])+f[1];
   u[n]=-h2i*(au[n-1]-2*au[n])+f[n];

   //输出残差的最大范数
   for (i=2;i<n;i++){
     mu[i]=u[i];
   }
   umaxnorm=maxnorm_one(a,mu,n);
   rateR=(pumaxnorm-umaxnorm)/pumaxnorm;
   std::cout<<"The maximum norm of residuals after "<<icycle<<" times FMG cycle="<<umaxnorm<<std::endl;
   std::cout<<""<<std::endl;
  std::cout<<"The "<<icycle<<"-th residual reduction rate is about "<<rateR<<std::endl;
  std::cout<<""<<std::endl;
  //判断精确度是否达到
   for (i=1;i<=n;i++){
   er[i]=y[i]-au[i];
   }
   Maxnorm=maxnorm_one(a,er,n);
  std::cout<<"The maximum norm of error after "<<icycle<<" times FMG cycle="<<Maxnorm<<std::endl;
  std::cout<<""<<std::endl;
  if (umaxnorm<RA) break;
 }
}
